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Strongly correlated systems with geometric frustrations can host the emergent phases of matter with uncon¬ 
ventional properties. Here, we study the spin S = 1 Heisenberg model on the honeycomb lattice with the 
antiferromagnetic first- (Ji) and second-neighbor (J 2 ) interactions (0.0 < J 2 /J 1 < 0.5) by means of density 
matrix renormalization group (DMRG). In the parameter regime J\ < 0.27, the system sustains a Neel anti¬ 
ferromagnetic phase. At the large J 2 side J 2 I J\ ^ 0.32, a stripe antiferromagnetic phase is found. Between the 
two magnetic ordered phases 0.27 < J 2 /Ji < 0.32, we find a non-magnetic intermediate region with a plaque- 
tte valence-bond order. Although our calculations are limited within 6 unit-cell width on cylinder, we present 
evidence that this plaquette state could be a strong candidate for this non-magnetic region in the thermodynamic 
limit. We also briefly discuss the nature of the quantum phase transitions in the system. We gain further insight 
of the non-magnetic phases in the spin-1 system by comparing its phase diagram with the spin-1/2 system. 


I. INTRODUCTION 


Since Anderson proposed the resonating valence bond the¬ 
ory to explain the high-temperature superconductivit}® the 
study of spin liquid (SLj^tS] in frustrated magnetic sy stems 
have been attracting much attentions for almost 40 year^SEH 
In recent years, this field has achieved exciting progresses 
with identifying realistic examples of SL states. Among the 
various SL candidates, the most promising candidate is the 
spin-1/2 kagome Antiferromagnet. In experimental side, the 
strong evidences supporting the gapless SL have been dis¬ 
covered in the^in-1/2 kagome antiferromagnet materials 
Herbertsmithit^^G^ and Kapellasit^i^l^. In theoretical stud¬ 
ies, the spin-1/2 kagome Heisenberg model has been found to 
sustain a SL ground state although the nature of the SL is still 
under debate between the gapped Z 2 SL®S] obtained from 
density-matrix renormalization group (DMRG) and the gap¬ 
less [/(I) Dirac SL favored in the variational studies of the 
Gutzwiller projected wavefunctiorP^®! Very recentl y, by in- 
troducing the second and third neighbor interactions^SEI or 
the chiral interaction^^ in kagome systems, a gapped chi¬ 
ral spin liquid which breaks time-reversal symmetry is un¬ 
ambiguously established as the 1 / = 1/2 fractional quan¬ 
tum Hall state through fully characterizing the topological 
propertieof the state. 

Besides the kagome systems, the spin-1/2 J 1 -J 2 Heisen¬ 
berg models on the square and honeycomb lattices have also 
being considered as the promising candidates of SL. In par¬ 
ticular, the SL on square lattice is considered significant to 
understand the high-temperature superconductivity in copper 
oxidd^l Recently, the long-debated non-magnetic regions in 
these two models have been studied intensively as the po ssi- 
ble realizations of gapped Z 2 SL^SElor gapless SlEMI By 
performing DMRG calculations on cylinder systemJ^^®, it 
is found that the Z 2 SL behaviors in the intermediate region 
of the both models appear not stable on the wide systems. In¬ 
stead, a plaquette valence-bond (PVB) state may dominate the 
non-magnetic regions. 

The frustrated spin-1 magnetic systems on the square and 
honeycomb lattices are also particularly interesting as they 
may be the parent magnetic systems for the iron-based su¬ 
perconductivity. The square Heisenberg models with frustrat- 
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FIG. 1: Quantum phase diagram of the spin-1 J 1 -J 2 Heisenberg 
model on the honeycomb lattice. With increasing J 2 coupling, the 
system has a Neel AFM phase for J 2 < 0.26, a stripe AFM phase for 
J 2 // 0.32, and an intermediate non-magnetic phase with a plaque¬ 
tte valence-bond order in our DMRG calculations. As the finite-size 
effects on different cylinder geometries, the first transition point is 
estimated at 0.26 ~ 0.28. 


ing further-neighbor couplings and quadratic interactionJ^^® 
have been studied intensively using the mean-field analysis to 
investigate the possible nematic order. By using DMRG cal¬ 
culations on the spin-1 J 1 -J 2 square Heisenberg modei^, a 
non-magnetic phase between the Neel and the stripe antifer¬ 
romagnetic (AFM) phase is obtained. However, the nature of 
this non-magnetic phase is far from clear although a nematic 
paramagnetic state has been proposed based on the field- 
theory descriptiorl^. On the other hand, the frustrated spin-1 
honeycomb Heisenberg models have not been studies system¬ 
atically, which may be relevant to the spin model for the hon¬ 
eycomb iron-based superconductivity m aterial Sr PtAp^^. 
For the spin-1/2 J 1 -J 2 honeycomb modeiSEUSIIllHll!^ a Neel 
AFM phase and a staggered dimer phase are found at J 2 < 
0.22 and J 2 > 0.35, respectively. Between these two phases 
for 0.25 < J 2 < 0.35, a PVB phase is identified in DMRG 
calculationPalUJas the PVB corre lation length keeps growing 
fast with system width on cylindei^ffi] por spin 5 = 1 J 1 -J 2 
model, the studies are rare and it is unclear whether the quan¬ 
tum phases such as the PVB and the staggered dimer phases 
would persist with spin magnitude increasing from 1/2 to 1, 
and what kinds of classical states might emerge in the honey¬ 
comb systerrl^. 

In this article, we study the spin-1 Heisenberg model on 
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(a) AC4-6 (b) ZC4-6 


FIG. 2: Cylinder geometries used in the DMRG calculations, (a) is 
an AC4-6 cylinder and (b) is a ZC4-6 cylinder. 


the honeycomb lattice with the frustrating J 1 -J 2 AFM in¬ 
teractions by using the DMRG with spin rotational SU{2) 
symmetrjESEil xhe Hamiltonian of the model is given as 

iF = Ji ^ S, • S, + J 2 ^ S, • Sj, (1) 

(ij> 

where Ji and J 2 are the first- and second-neighbor AFM inter¬ 
actions. We set Ji as the energy scale, and lattice spacing be¬ 
tween the nearest-neighbor sites as the length scale. Through 
our SU (2) DMRG calculations on cylinder systems, we es¬ 
tablish a quantum phase diagram as shown in Fig. B By 
studying the spin correlation function, we find a Neel AFM 
phase for 0 < J 2 < 0.27. For 0.5 > J 2 ^ 0.32, we find the 
magnetic ordered state rather than the staggered dimer state. 
In this region, the obtained magnetic order state depends on 
the cylinder geometry in our finite-size calculations. By com¬ 
paring the ground-state bulk energy on different cylinder ge¬ 
ometries, we find that the stripe AFM state always possesses 
the lower energy and thus appears to be the ground state in 
the thermodynamic limit. Between the two magnetic order 
phases 0.27 J 2 ^ 0.32, the system has a narrow non¬ 
magnetic region with the non-uniform bond energy on wide 
cylinders. During the increase of the kept states from 2000 
to 8000 SU(2) states, a PVB dimer order is found stabilized 
on the studied cylinder systems, which suggests that the PVB 
state is a strong candidate for this intermediate phase region. 
Finally, we discuss the nature of the quantum phase transi¬ 
tions in the system with the help of the bipartite entanglement 
entropy. 

In our DMRG calculations, we study the cylinder systems 
with width up to 8 (6) unit cells in the magnetic ordered phases 
(intermediate phase), by keeping up to 8000 5(7(2) states to 
ensure the convergence. The truncation error is controlled be¬ 
low 10“® for Ly = A (Ly is the number of unit cell in the y 
direction) cylinder and below 10“^ for the other calculations. 
The cylinder geometries are shown in Fig.|^ The first cylinder 
ACm-n has the armchair open edges, where m is the number 
of two-site unit cells along the y direction and n is the number 
of columns along the x direction. The second cylinder has the 
zigzag open edges and is denoted as ZCm-n cylinder. 


(a) J2=0.0, AC4-24, spin correlations 



(b) J2=0.0, ZC4-24, spin correlations 



FIG. 3: Spin correlation functions in real space for J 2 = 0.0 on 
(a) AC4-24 cylinder and (b) ZC4-24 cylinder. Tbe green site is tbe 
reference site in tbe middle of cylinder, and tbe blue solid and red 
shaded circles denote tbe positive and negative correlations, respec¬ 
tively. Tbe radius of circle is proportional to tbe magnitude of corre¬ 
lations. Here, we only show tbe left balf 2 x 4 x 12 sites. 

II. NEEL AEM PHASE 

First of all, we study the Neel AFM phase in the small J 2 
side. Due to the limit of system width in the DMRG calcu¬ 
lations for the spin-1 system, we do not have enough data of 
magnetization on different system widths to estimate the re¬ 
sult in thermodynamic limit through extrapolation. Instead, 
on the finite-width cylinders, we calculate the spin correlation 
functions along the cylinder axis direction (the x direction) 
and study their decay behaviors with increasing J 2 . 

In Fig. we demonstrate the spin correlation functions in 
real space for J 2 = 0.0 on the AC4-24 and ZC4-24 cylin¬ 
ders. Clearly, the spin correlations exhibit a Neel AFM pat¬ 
tern with two magnetic sublattices. We follow the J 2 depen¬ 
dence of the spin correlation decaying to detect the vanishing 
of Neel order. On the AC4-24 cylinder in Fig. Qa), we find 
that the spin correlation length keeps decreasing with grow¬ 
ing J 2 , which reaches a minimum at J 2 ~ 0.25. Slightly 
above J 2 = 0.25 such as J 2 = 0.27 as shown in Fig.j^a), the 
Neel AFM pattern of correlations is destructed, which signals 
a phase transition with vanishing the Neel order. On the ZC4- 
24 cylinder as shown in Fig. Qb), we find that the spin cor¬ 
relations decay much slower than those on the AC4-24 cylin¬ 
der near the transition point, which indicates the finite-size 
effects of the system when approching phase boundary. Be¬ 
yond J 2 ^ 0.28, the spin correlations decay fast as shown in 
Fig.j^b) at J 2 = 0.29. Based on the spin correlations on both 
AC4 and ZC4 cylinders, we estimate the Neel order vanishing 
at J 2 ~ 0.26 - 0.28. 

Based on the decay behaviors of spin correlations, we 
plot the J 2 dependence of the long-distance spin correlations 
Sd = \/|(5o5£;)| {d is the longest distance in Figs, [^a) 
andj^b)) as shown in Figs.Qc) andQd). We find that the 
vanishing of Sd and the destruction of the Neel AFM patter 
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FIG. 4: (a) and (b) are the log-linear plots of spin correlations on the 
AC4-24 and ZC4-24 cylinders, (c) and (d) are the J 2 dependence of 
the long-distance spin correlations Sd = \/|(5'oS'd)| (d is the longest 
distance in (a) and (b)). In the subfigure (a), the correlation data for 
J 2 < 0.27 are all positive. For J 2 = 0.27, the Neel AFM pattern is 
destructed, thus the correlations change sign in some places. 


on both geometries are consistent with a phase transition at 
J 2 ~ 0.26 — 0.28. We have also checked the spin correlations 
on AC6 and ZC6 cylinders. Although the long-range corre¬ 
lations on these wider cylinders are not fully converged, their 
behaviors are qualitatively consistent with those on the AC4 
and ZC4 cylinders. 


III. STRIPE AFM PHASE 

In the spin-1/2 J 1 -J 2 honeycomb Heisenberg model, the 
system is in a staggered dimer phase with breaking lattice ro¬ 
tational symmetry and short-range spin correlations for J 2 > 
0.35. Interestingly, in this spin-1 system we hnd magneti¬ 
cally ordered states instead of the staggered dimer. As shown 
in Fig. I^a) for J 2 = 0.4 on the AC6-18 cylinder, we find 
a magnetic ordered state with the 8-site unit cell denoted by 
the green dashed rectangle. This state is stable on all the AC 
cylinders that we have studied (AC4, AC6, and ACS). How¬ 
ever, on the ZC6-18 cylinder at J 2 = 0.4, we find a stripe 
AFM state as shown in Fig. |^b). This stripe state is also sta¬ 
ble on the different ZC cylinders (ZC4, ZC6, and ZC8). In 
Figs. a) andj^b) of the log-linear plots of spin correlation 
functions at the large J 2 side, we find that the spin correlation 
length diminishes with decreasing J 2 . On AC4 cylinder, the 


(a) J2=0.4, AC6-18, left half lattice 



(b) J2=0.4, ZC6-18, left half lattice 



FIG. 5: Spin correlation function {So ■ Sj) in real space for J 2 = 0.4 
on (a) AC6-18 cylinder and (b) ZC6-18 cylinder. The green site is 
the reference spin So in the middle of lattice, and the blue solid and 
red shaded circle denote the positive and negative spin correlations, 
respectively. The area of the circle is proportional to the amplitude 
of correlations. The dashed rectangles denote the unit cells. 

spin correlations have a sharp increase for J 2 > 0.35; and on 
ZC4 cylinder, correlations grow rapidly for J 2 > 0.32. For 
J 2 < 0.32, the spin correlations decay quite fast to vanish, 
which is consistent with a non-magnetic phase region. 

To identify which state is the exact ground state at large 
J 2 side, we compare the ground-state energy on both cylin¬ 
der geometries as shown in Fig. |^c) for J 2 = 0.4. We ex¬ 
tract the ground-state energy from the bulk bond energy on 
long cylinder systems. Interestingly, we find that the energies 
on ZC cylinder are always lower than those on AC cylinder. 
As we expect, the two geometries should give the same en¬ 
ergy in the thermodynamic limit. Thus, the different energy 
in Fig. I^c) indicates the strong hnite-size effects on systems 
we can study. Although we cannot definitely determine which 
geometry gives the correct ground state in the thermodynamic 
limit, based on our calculations we believe that the stripe state 
on ZC cylinder would win because it always has the lower 
energy for system sizes we studied. 

As shown in Fig. of the spin configuration in the stripe 
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FIG. 6: Log-linear plots of spin correlations at large J 2 side for (a) 
AC4-24 and (b) ZC4-24 cylinders, (c) Cylinder width dependence of 
the bulk ground-state energy for J 2 = 0.4 on the AC and ZC cylin¬ 
ders. ZC cylinders always have the lower energy than AC cylinders. 


J2=0.5, ZC4-24, bond energy 

0.699 0.748 0.752 



FIG. 7: The nearest-neighbor bond energy for J 2 = 0.5 on the ZC4- 
24 cylinder. Here, we only show the left half lattice. The positive 
horizontal (labeled by blue) and the negative vertical (labeled by red) 
bond energies are consistent with the spin configuration in the stripe 
AFM state. 


AFM phase, we expect that the nearest-neighbor bond energy 
{Si ■ Sj) of the two parallel spins are positive and those be¬ 
tween the anti-parallel spins are negative. This feature also 
can be used to characterize the appearance of the stripe phase. 
In Fig. 1^ we demonstrate the bond energy for J 2 = 0.5 on 
ZC4-24 cylinder, which indeed has the bond energy pattern 
expected for the stripe AFM phase. With decreasing J 2 , we 
find that this bond energy pattern could persist to J 2 = 0.32, 
which is consistent with the phase transition point estimated 


from spin correlation function. 


IV. INTERMEDIATE PHASE REGION 


Next, we investigate the intermediate region between the 
two magnetic ordered phases. First of all, we calculate the 
first-neighbor bond energy (Si-Sj) to study the possible lattice 
symmetry breaking. To accommodate the possible valence- 
bond solid states on cylinder system, we perform calcula¬ 
tions on two cylinder geometries, the AC cylinder and the 
trimmed ZC (tZC) cylinder with some tri mmed sites on the 
open boundaries (see the lattice in Fig.[^b)j^2lliJ, On the AC4 
and tZC4 cylinders, the bond energies are quite uniform in the 
bulk of cylinder, which is consistent with the short order cor¬ 
relation length on narrow system. Thus, to detect the possible 
lattice symmetry breaking, we need go to the wider systems to 
study the behavior of correlation length with growing cylinder 
widtlPSEl. As the DMRG convergence in the intermediate re¬ 
gion is quite challenging on the wider systems, we can only 
study the AC6 and tZC6 cylinders. During the DMRG calcu¬ 
lations, we measure the bond energy by increasing the optimal 
state number step by step. By keeping the states up to 8000 
5(7(2) states, we find a PVB dimer pattern stabilized as shown 
in Fig.[^ which is robust in the whole lattice and thus has bro¬ 
ken the lattice translation symmetry on the finite-size cylin¬ 
ders. For the calculations on the AC6 and tZC6 cylinders by 
keeping 8000 states, the trTincation errors are about 1 x 10“®. 
Although we cannot obtain the more accurate results, based 
on the clearer PVB pattern with growing kept state number, 
we argue that the PVB order is stable on these geometries. In¬ 
terestingly, the field-theory analyses have indicated that such 
a PVB state could emerge proximate to the Neel AFM phase 
on both the spin-1/2 and the spin-1 honeycomb Heisenberg 
model^. We also study the bond energy for other J 2 around 
the intermediate region. For J 2 = 0.31, the systems also ex¬ 
hibit the PVB pattern. However, for J 2 > 0.32 the bond en¬ 
ergy pattern appears like Fig. in consistent with the stripe 
state. 

To further distinguish the three different phases, we cal¬ 
culate the spin gap on long cylinder. As the broken 5(7(2) 
symmetry and emerging Goldstone boson, the spin gap is ex¬ 
pected to vanish in both the Neel AFM and the stripe AFM 
phases. However, as the spin singlet bonds formed in the 
valence-bond solid states, the spin gap should be non-zero 
in a PVB state. To find the spin gap on cylinder, we follow 
the standard method in DMRG calculations. We first sweep a 
long cylinder to find the ground state with the energy Eq , and 
then we sweep the bulk by targeting the total spin 5 = 1 sec¬ 
tor to find the lowest-energy state with the energy Ei. Then 
the spin gap of the bulk At is the difference between Eq and 
El, At = El — Eq. We calculate the long cylinder to make 
sure the spin gap At is converged and independent of cylin¬ 
der length. We do not show the data on AC cylinders as their 
ground state in the large J 2 side has the higher energy on our 
studied systems. As shown in Fig. |^a), the spin gap At is 
vanishing-small for < 0.25, which is consistent with the 
Neel AFM order. At J 2 ~ 0.27, the spin gap increases dra- 






5 


(a) J 2 = 0.3, AC6-24, M = 8000 



(b) J2 = 0.3, tZC6-18, M = 8000 

- 0.051 - 0.066 - 0.069 - 0.069 - 0.069 



FIG. 8: The nearest-neighbor bond energy textures for J 2 = 0.3 on (a) AC6-24 cylinder and (b) tZC6-18 cylinder, which are obtained by 
keeping 8000 SU (2) states. The bond textures are obtained by subtracting the average bond energy. The blue and red bonds denote the positive 
and negative bond textures, respectively. On both geometries, a PVB bond pattern is found. For AC6-24 in (a), only the left half lattice is 
shown here. 


matically and then decreases for J 2 > 0.35. The enhanced 
spin gap in the intermediate region could be served as an ev¬ 
idence for the valence-bond solid states. In the stripe phase 
J 2 ^ 0.32, the spin gap is large on ZC4 cylinder, which is re¬ 
lated to finite-size effects because it should be vanished in the 
thermodynamic limit. This size effect can be understood from 
the bond energy in Fig.|^ which has the AFM bonds along the 
y axis. Thus, the spin ttiplet gap is strongly correlated with the 
size scale in the y direction, which accounts for the large gap 
on the narrow ZC4 cylinder. To show the size dependence of 
spin gap on the cylinder width, we also calculate the spin gap 
on the ZC6 cylinder deep inside the stripe phase. As expected, 
the spin gap on ZC6 cylinder drops rapidly, consistent with the 
vanished gap in the thermodynamic limit. We do not show the 
spin gap on ZC6 cylinder for the intermediate region because 
the DMRG calculations in the spin-1 sector are quite far from 
convergence in this region. 

Finally, we demonstrate the bulk ground-state energy for 
J 2 = 0.3 on different cylinders. We obtain the bulk energy 
from calculating all the bond energy in the middle of cylinder. 
As shown in Fig. [9jb), the ground-state energy at J 2 = 0.3 
smoothly increases with growing cylinder width, which be¬ 
haves differently from the energy scaling in the stripe phase 
demonstrated in Fig. |^c), where the energies on the two ge- 



FIG. 9: (a) J 2 dependence of the spin gap on the ZC4 and ZC6 cylin¬ 
der systems, (b) Cylinder width dependence of the bulk ground-state 
energy on AC (AC4 and AC6) and ZC (ZC4 and ZC6) cylinders for 
J 2 = 0.3. 


ometries scale separately. We also notice that the ground-state 
energy per site at J 2 = 0.3 changes slightly with growing sys¬ 
tem width. The energy appears to approach Coo — —1.262, 
which provides an upper bound for the ground state energy in 
the thermodynamic limit. 
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FIG. 10: (a) J 2 dependence of the bipartite entanglement entropy on 
the ZC4 cylinder. At J 2 ~ 0.28 and 0.32, the entropy has the sharp 
jump and drop, which appear consistent with the quantum phase tran¬ 
sitions. Entanglement spectra for (b) J 2 = 0.0 on ZC6-18 cylinder, 
(c) J 2 = 0.4 on ZC6-18 cylinder, and (d) J 2 = 0.3 on ZC4-24 cylin¬ 
der. The blue numbers denote the number of the largest eigenvalues 
in each S sector. The green dashed lines denote the tower of states 
structure in the magnetic ordered states. 


V. ENTANGLEMENT ENTROPY AND SPECTRUM 


To further characterize the different phases and phase tran¬ 
sitions, we study the bipartite entanglement entropy and en¬ 
tanglement spectrum. In Fig. [T0|a), we demonstrate the J 2 
dependence of entropy on the ZC4 cylinder. The entropy 
shows a sharp jump at J 2 — 0.28 and a drop at J 2 ~ 0.32, 
which have been shown to characterize various phase transi¬ 
tions in one-dimensional systemJ^HEll and are consistent with 
the identified transitions points. While the discontinuous jump 
at J 2 ~ 0.28 suggests a first-order transition from the Neel 
AFM to the non-magnetic phase, the smoother entropy de¬ 
crease near J 2 — 0.32 might be consistent with a weak hrst- 
order transition. 


In Figs. 10 b-d), we demonstrate the bipartite entanglement 


spectrum in each phase. For the ordered phases with con¬ 
tinuous symmetry breaking, the lower part of the entangle¬ 
ment spectrum is in correspondence with the “tower of states” 
(TOS) spectrunJ221E21, in Figs. 10 h) and 10 c) of the spec¬ 
tra in both the magnetic ordered states with breaking SU{2) 
to U{1) symmetry, the spectra have a single dominant eigen¬ 
value in each S sector, which are separated from the higher 


levels by the entanglement gap and follow a linear behavior 
with S{S + 1) (S is the quantum number of total spin). All 
these features are consistent with the TOS structures of the 
energy spectra in the corresponding magnetic ordered states 
on the honeycomb latticd^. In the intermediate phase region, 
the spectrum is totally different from the magnetic ordered 
states, as shown in Fig.fl^d) for J 2 = 0.3 on ZC4 cylindeJ^. 
The low-lying degeneracy in the S' = 1 sector changes from 
1 to 2, and there is no clear entanglement gap between the 
largest eigenvalues and the rest of spectrum. These features 
distinguish the intermediate phase from the magnetic ordered 
phases. 


VI. SUMMARY 


We have studied the quantum phase diagram of the spin- 
1 J 1 -J 2 Heisenberg model on the honeycomb lattice using 
density-matrix renormalization group calculations on cylinder 
system. We have established three different phases including 
two magnetic ordered phases and a non-magnetic phase. For 
J 2 < 0.27, we find a Neel AFM phase. For J 2 > 0.32, we 
find two possible candidate magnetic ordered states depend¬ 
ing on the different geometries. On AC cylinders, the system 
is a magnetic order state with the 8-site unit cell, while on 
ZC cylinders it is a stripe AFM state. By comparing the bulk 
ground-state energy on the two geometries, we find that the 
ZC cylinders always have the lower energy than the AC cylin¬ 
ders, which strongly suggests the stripe AFM state as the true 
ground state in the thermodynamic limit. 

Between these two magnetically ordered phases with 
0.27 ^ J 2 ^ 0.32, we hnd a non-magnetic phase region. 
On both AC6 and tZC6 cylinders, the systems have the 
non-uniform bond energy. By increasing the kept states to 
8000 SU{2) states (equivalent to about 24000 U{1) states), 
we find a stable plaquette valence-bond order emerging in the 
systems. The spin gap on finite-size cylinder also enhances 
dramatically in this phase region. Our results indicate that 
the plaquette state is a strong candidate for this non-magnetic 
phase. Moreover, the sizable entropy change on the phase 
boundaries indicates that the nature of the phase transitions 
from the magnetic ordered phases to the non-magnetic phase 
might be first order. 
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